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The cluster mean-field approximations are performed, up to 13 cluster sizes, to study the critical 
behavior of the driven pair contact process with diffusion (DPCPD) and its precedent, the PCPD 
in one dimension. Critical points are estimated by extrapolating our data to the infinite cluster size 
limit, which are in good accordance with recent simulation results. Within the cluster mean-field 
approximation scheme, the PCPD and the DPCPD share the same mean-field critical behavior. The 
application of the coherent anomaly method, however, shows that the two models develop different 
coherent anomalies, which lead to different true critical scaling. The values of the critical exponents 
for the particle density, the pair density, the correlation length, and the relaxation time are fairly 
well estimated for the DPCPD. These results support and complement our recent simulation results 
for the DPCPD. 
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The absorbing phase transition (APT) has been stud- 
ied extensively to understand many-body cooperative 
phenomena in nonequilibrium systems Up to now, 
two universality classes have been firmly established: di- 
rected percolation (DP) and parity conservation (PC) 
universality classes 2]. A few other candidates for dif- 
ferent universality classes have been reported in recent 
literatures. One is the DP system coupled with a static 
conserved field 0]. Although the reported values for the 
critical indices are rather scattered 'sl, , it is widely 
believed that these systems form a universality class, dif- 
ferent from the DP and the PC class. Another candidate 
is the pair contact process with diffusion (PCPD) that 
has as yet defied any consensus on the universality issue. 
Various scenarios have beenproposed, including a new 
single universality class 0, 13, a marginally perturbed 
DP process with continuously varying exponents 0| , and 
a DP process with a huge crossover time 0, ^| , which 
are summarized in a recent review |l2j. 

Recently, we studied the driven PCPD (DPCPD) 
which is a variant of the PCPD by introducing biased 
diffusion It is shown that the driving is relevant 

and the DPCPD exhibits a "mean-field-like" critical be- 
havior even in one dimension. Since the DP class is in- 
sensitive to the driving, the DP scenario with a huge 
crossover time should be eliminated. There was a recent 
attempt to understand the PCPD using the renormaliza- 
tion group (RG) analysis on a single-species Bosonic ac- 
tion derived from the microscopic master equation. How- 
ever, it turned out to be improper to describe the critical 
behavior of the PCPD __14J. In our previous work [isj . 
we pointed out a possible reason for this failure and sug- 
gested that the PCPD may be described properly by a 
field theory with two independent fields. Still, the search 
for the coarse-grained action adequate for the PCPD re- 
mains a challenge. 

Besides the RG technique on the proper action |[l5l |. 
there are a few other efficient methods to investigate 
the absorbing critical phenomena. Numerical simulations 



along with a finite-size-scaling analysis [ig and direct in- 
tegrations of corresponding Langevin equations [l7j | are 
two typical examples. Another frequently used method 
is the cluster mean- field (CMF) approximation [l3| fol- 
lowed by the coherent-anomaly method (CAM) analysis 
[l^ . This method is known to be effective to obtain a 
quantitative phase diagram and sometimes even explore 
a true critical scaling behavior However, the ac- 

curate measurement of critical indices is only limited to 
rather simple DP systems. More complex critical behav- 
iors like in the PC and the PCPD classes could not have 
been probed with a reasonable accuracy as yet by the 
CAM analysis IMIll- 

In this paper, motivated by our recent results that the 
DPCPD exhibits a distinct critical behavior from the 
PCPD and also a mean-field-like behavior even in one 
dimension ^13] , we develop the CMF approximations for 
the DPCPD and the PCPD, expecting that the CAM 
analysis would produce a reasonable estimate for the 
mean-field-like critical indices of the DPCPD. Also, di- 
rect comparison of the CMF data for two models may 
provide an independent support for different scaling be- 
haviors. 

We set up dynamic CMF equations up to ti = 13 
cluster size. The steady-state solutions are obtained 
within machine accuracy using MATHEMATICA. Dynamic 
information is also extracted from the smallest eigen- 
value of the "linearized" transition matrix. Subsequently, 
through the CAM analysis, we estimate the values of the 
critical exponents for the particle density, the pair den- 
sity, the correlation length, and the relaxation time. 

The model is defined on a one-dimensional lattice of 
L sites with periodic boundary conditions. At each site, 
there is at most one particle and no multiple occupancy 
is allowed. Hence the configuration is specified by the 
occupation number which is either 1 or at every site. 
Each particle hops to the right (left) with transition rate 
Dr {Dl)- The total number of particles in the system 
varies by branching and annihilating events mediated by 
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a particle pair {2 A 3 A and 2 A — » 0). The transition 
rate is p (1 — p) for the annihilating (branching) event 
with < p < 1. These three dynamics can be described 
by the master equation which takes the form 



10-' 



-H\P;t), 



(1) 



where \P; t) is the state vector at time t and the "Hamil- 
tonian" is written as H = X)i=i with 

Hi = DR{piVi+i ~ CLidl^^) + DL{viPi+i - a,|a,+i) 
- p{diai+i - pipi+i ) (^2) 

where ai(a|) is the annihilation (creation) operator of 
hard core particles, satisfying {0^,0^} — 1 and [0.^,0^] — 
[aijOj] — for i ^ j, pi ~ alai is the number operator, 
and Vi = 1 — Pi. 

Three different cases arise depending on the values of 
Dji and Dl- The case oi Dn = Dl = represents the 
pair contact process (PCP) which has infinitely many ab- 
sorbing states and is known to belong to the DP class at 
least for static situations [23 ■ Since the cluster approxi- 
mations along with the CAM analysis have been already 
performed for the PCP by several authors previously |24l| . 
we skip the analysis of the PCP here. The PCPD corre- 
sponds to Da = Dl ^ and the DPCPD corresponds 
to Da ^ Dl. In what follows, we set Da + Dl — 1 
for convenience and Dn is chosen to be 1/2 (1) for the 
PCPD (DPCPD). 

We consider an n-site probability function P„(p;t). It 
is defined as the probability at time t to find an n-site 
cluster of the configurational state p = (pi,p2, ■ • ■ ,Pn), 
where an occupational state pi at site i takes either or 
1. Tracing out Eq. over occupational states outside 
the cluster {{pi} with i < or i > n -|- 1), one may find 
a formal exact expression 



dP„{p;t) 
dt 



Fp{Pn, Pn+l, Pn+2) , 



(3) 



where the function Fp involves the sets of n-, (n -I- 1)-, 
and [n + 2)-site probability. Notice that Pn+i and Pn+2 
terms show up due to the boundary dynamics of the n- 
site cluster. 

As the infinite hierarchy appearing in Eq. Q is the 
major obstacle towards analytic treatment, we need an 
approximation scheme to truncate the hierarchy at finite 
n. In this paper, we take the so-called (n -I- \,n) ap- 
proximations [13, where P„+2 and Pn+i are expressed 
in terms of products of P„'s. Then, the rate equations 
for the n-site cluster probability function become 



dPnjp] t) 
dt 



= Fp{{Pn}). 



(4) 



where P„ is now the approximate (mean-field) probabil- 
ity function. 
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FIG. 1: Log-log plot of pB vs p^^ — p obtained from the 12- 
cluster CMF approximation for the DPCPD. The slope of the 
straight line is 1. In the inset, pp's are plotted against pj'^ —p. 
The slope of the straight line is 2. 



The stationary probability distribution function P^(p) 
can be obtained by solving the set of coupled equa- 
tions Fp = 0. For given n, the number of equations 
and the number of variables are both 2", but not all 
are independent. The translational invariance and the 
normalization condition guarantee that all Pn{p) with 
pi = can be expressed in a linear combination of 
P„(p)'s with pi = 1. For example, in case of n = 4, one 
can easily show that P4(0011) P3(011) - P4(1011) = 

-P2(ll) - P3(lll) - P4(1011) - P3(110) - P4(1011) - 

P4(1101) + P4(1100) - P4(1011). The absorbing (vac- 
uum) probability P4(0000) can be determined by the nor- 
malization condition. Hence the DPCPD case has 2"^^ 
independent variables. In case of the PCPD, the left- 
right symmetry further reduces the number of indepen- 
dent variables, for example, P4(1101) = P4(1011) and so 
on. 

We use MATHEMATICA to find the stationary solutions 
Pnip) up to n = 13 with machine accuracy (10~^°) for 
given values of parameters, p and Dfi With P*, we 
calculate the particle density ps and the pair density pp 
in the steady state as 

Ps - P^(l) = TrpP^(p)<5p,,i, 

Pp = P|(ll) = Trp P,:(p)(5,„i<5p,^,,i, (5) 

where k (and fc -I- 1) denotes an arbitrary site inside the 
cluster. At a fixed value of Djj, the order parameters, 
Ps and Pp, simultaneously vanish for large p (pair annihi- 
lation rate) and the system exhibits an absorbing phase 
transition into vacuum at p = p". 

Near the transition point p", the order parameters 
scale as 

^MF /^MF 

Ps^A,M-pf' , Pp ^ B^{p: - p)"^ , (6) 

where we find the mean-field values for the order param- 
eter exponents: = 1 and (3^^ = 2. Figure ^ shows 



3 



TABLE L Cluster approximation results for the PCPD 
model. The errors are in the last digits. 
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the n = 12 cluster approximation results for the DPCPD. 
We estimate the critical point p" and the critical ampli- 
tudes An and i?„ by fitting five data near the transition 
{\p^p2\ < 5 X 10~^), linearly for ps and quadratically for 
Pp. Our results are tabulated in Table U for the PCPD 
model and in Table|Tllfor the DPCPD model. Notice that 
the relative errors for are extremely small (~ 10~^), 
but the amplitudes A„ and _B„ still have a sizable relative 
error (~ 10^"*). 

It is interesting to note that, for n < 3, the diffusion 
bias does not enter the CMF rate equations at all. The 
functional form of Fp{{Pn}) in Eq. ^ is identical for the 
PCPD and the DPCPD. The left-right symmetry among 
Pnipys, is automatically enforced due to the translational 
invariance, regardless of the details of the dynamics. For 
example, ^3(110) = ^2(11) - ^3(111) = ^3(011) and so 
on. However, for n > 4, the translational invariance does 
not guarantee the left-right symmetry, which may be bro- 
ken by the dynamics with a broken left-right symmetry. 

In nonequilibrium systems, dynamic relaxation behav- 
ior provides one of the key pieces of information on the 
system. Off criticality, the order parameters are expected 
to approach their stationary values exponentially with a 
characteristic relaxation time r. At criticality, r diverges 
and the order parameters decay algebraically. One may 
roughly estimate r by numerically integrating the rate 
equations (@J and fitting time-dependent data into an 



TABLE II: Cluster approximation results for the DPCPD 
model. The errors are in the last digits. 
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n 
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0.216 140 3513 
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44.37 
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0.194 381 7410 
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0.180 689 8311 
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10 


0.177 954 3360 


4.410 


103.7 
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11 


0.175 711 7674 


4.397 


112.3 


17.97 


12 


0.173 837 8803 


4.383 


120.8 


17.99 


13 


0.172 247 8976 


4.368 


128.9 





FIG. 2: Log-log plots of A„ vs 1/n for the PCPD and 
the DPCPD. In this figure, the value of for the PCPD 
(DPCPD) is set to be 0.1335 (0.153). 



exponential form. However, this method does not pro- 
duce high-precision data. In this paper, we propose a 
different method to calculate r with machine accuracy in 
the CMF approximation scheme. 

Since the stationary solutions of Eq. Q were obtained 
with machine accuracy, we can linearize Eq. I^J near the 
stationary solutions very accurately. The linearized equa- 
tion takes the form 



d\Pn;t) 

dt 



= -M\Pn;t), 



(7) 



where |P„; i) is the (n-cluster) state vector with the com- 
ponents Pn (p; t) and M is a square matrix. It is trivial to 
show that the eigenvalues of M are equal to the inverse 
of various characteristic time scales of the dynamics. The 
most dominant slow mode is determined by the smallest 
eigenvalue As, i.e., the relaxation time r = Aj^. 

We analyze the linearized equation up to n = 12. Near 
criticality, we find 



(8) 



where we find again the mean-field value for the relax- 
ation exponent, i^^^ = 2. We estimate the critical points 
independently, which are found to be consistent with 
previous estimates from the density data in Tables and 
im where we also tabulate the estimated values for the 
ampHtude C„ for both the PCPD and the DPCPD. 

Now, we employ the coherent-anomaly method (CAM) 
introduced by Suzuki and co-workers [l9j and estimate 
the values of the true critical exponents. Following the 
CAM analysis, the n dependence of the critical point p" 
is predicted in the large n limit as 



A" 



(9) 



where A„ — p" — pc is the distance of p" from the true 
critical point Pc = lim„^oo Pc ^^'^ ^i- is the true correla- 
tion length exponent. 
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FIG. 3: The CAM analysis for the order parameters for the 
PCPD modeL The slope of the straight line is —1.51, which 
leads to (52 ~ 0.49. In the inset, the slope of the straight line 
is -0.084, which leads to /3i ^ 0.92. 



FIG. 4: The CAM analysis for the order parameters for the 
DPCPD model. The slope of the straight line is —0.92, which 
leads to 132 « 1.08. In the inset. An vs An is drawn without 
a log scale. An remains nearly constant which implies I3\ ~ 

= 1. 



We estimate Pc by applying the Bulirsch and Stoer 
(or BST) algorithm [^j to the series of {p"} and find 
that Pc = 0.134(2) for the PCPD and Pc = 0.154(3) for 
the DPCPD, which are in good agreement with simula- 
tion results of 0.133 522(2) and 0.151 032(1) [3. Al- 
ternatively, we estimate Pc and simultaneously using 
Eq. In Fig. [S] we plot A„ vs 1/n in a log- log plot, 
varying pc to find the best power-law fit. For the PCPD, 
the choice of Pc ~ 0.1335 yields the smallest fitting er- 
ror with v± = 1.04, where the data from n = 8 to 13 
are used. For the DPCPD, the best choice is Pc — 0.153 
with v±_ = 1.01. The relative error for pc is ~ 2%, and 
the error for v±_ is ~ 10%. The best simulation result of 
j/j. = 1.09(2) for the PCPD is within the errors, but the 
accurate measurement seems to be out of reach with data 
up to n = 13. Our estimate of v± = 1.01 for the DPCPD 
is in very good agreement with the expected mean-field 
value vf^ ^l. 

The amplitudes An and i?„ are expected to scale as 

Anr^ ^n^^"''-^'\ Bn^A-^^^"-^'\ (10) 

where Pi and P2 are the true critical exponents for the 
order parameters. In Fig. 13 the log-log plots of An and 
Bn vs An for the PCPD are presented. Here we use 
Pc ~ 0.133 522 (the best estimate from Monte Carlo sim- 
ulations) T3||. The CAM analysis leads to /3i « 0.92 and 
P2 ~ 0.49 [23, both of which are far from the simula- 
tion results of (ii k, (i2 ~ 0.36(2) 8, 9]. In particular, 
there is a huge discrepancy between the estimated values 
of j3i and /32 by the CAM analysis, which warns us that 
the CAM estimates for the order parameter exponents 
should be interpreted with great caution. This huge dis- 
crepancy also implies that the cluster sizes up to n = 13 
are still too small for the PCPD to reach the asymptotic 
regime where the system is dominated by long spatial cor- 
relations, induced by the long-term memory mediated by 
solitary particles |^. 



On the other hand, the CAM analysis for the DPCPD 
looks consistent with the simulation results. In Fig. 01 
we use Pc — 0.151 032 First, An seems not diverg- 
ing as A„ — > and reaching a nonzero constant, which 
implies fi\ — — 1. Second, Bn behaves very differ- 
ently from An and diverges with the exponent 0.92, 
which imphes that^= - 0.92 « 1.08. Numerical 
simulation results [l3j are in complete agreement with 
our CAM results. One should notice that (i^ does not 
assume the MF value, but seems to be equal to (5\ ex- 
cept a probable multiplicative logarithmic correction as 
found in the exponent iilv\\ by numerical simulations 
p^. This mean- field- like behavior is expected for the 
two-dimensional PCPD [28] , of which the upper critical 
dimension is believed to be 2. Our CAM results indepen- 
dently support the conclusion drawn from our numerical 
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FIG. 5: The CAM analysis for the relaxation time for the 
PCPD and the DPCPD. The slope of the straight line for the 
PCPD (DPCP) is -0.2 (-0.02) which leads to i/y = 1.8 (1.98). 
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simulations results [13 that the DPCPD critical behav- 
ior is distinct from the PCPD behavior and the upper 
critical dimension for the DPCPD is 1 rather than 2. 
Finally, we estimate from the relation 



-("11 



(11) 



where i^y is the true critical exponent for the relaxation 
time. In Fig. the log-log plots of C„ vs A„ are shown. 
We estimate that « 1.8 for the PCPD and i^^ « 1.98 
for the DPCPD. Rather surprisingly, the PCPD result 
is consistent with the simulation result of = 1.85(10) 
0]. For the DPCPD, the value of is quite close 



to the mean-field value of ly^^ 
simulation results. 



2, consistent with the 



In summary, we estimated the critical exponents for 
the PCPD and the DPCPD, using CMF approximations 
along with the CAM analysis. For the PCPD, the val- 
ues of the order parameter exponents are poorly esti- 
mated, while the estimates for the correlation and the 
relaxation exponents are consistent with simulation re- 
sults within error bars. In contrast, the CAM estimates 
for the DPCPD are in excellent accord with simulation 
results, supporting our conjecture that the upper critical 
dimension of the DPCPD is 1. 
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